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Abstract 

In the study of open quantum systems, one typically obtains the decoherence dynamics by solving 
a master equation. The master equation is derived using knowledge of some basic properties of 
the system, the environment and their interaction: one basically needs to know the operators 
through which the system couples to the environment and the spectral density of the environment. 
For a large system, it could become prohibitively difficult to even write down the appropriate 
master equation, let alone solve it on a classical computer. In this paper, we present a quantum 
algorithm for simulating the dynamics of an open quantum system. On a quantum computer, the 
environment can be simulated using ancilla qubits with properly chosen single-qubit frequencies 
and with properly designed coupling to the system qubits. The parameters used in the simulation 
are easily derived from the parameters of the system+environment Hamiltonian. The algorithm 
is designed to simulate Markovian dynamics, but it can also be used to simulate non-Mar kovian 
dynamics provided that this dynamics can be obtained by embedding the system of interest into 
a larger system that obeys Markovian dynamics. We estimate the resource requirements for the 
algorithm. In particular, we show that for sufficiently slow decoherence a single ancilla qubit could 
be sufficient to represent the entire environment, in principle. 

PACS numbers: 03.67.Ac, 03.67.Lx 
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I. INTRODUCTION 



For most problems of practical interest, the Schrodinger equation describing the evolution 
of a quantum system is too complicated to be solved exactly. Numerical simulations are 
therefore commonly employed for this purpose. However, simulating quantum systems on 
a classical computer is a hard problem. The dimension of the Hilbert space of the system 
increases exponentially with the size of the system (e.g., the number of particles in the 
system). On a quantum computer, on the other hand, the number of qubits required to 
simulate the system increases linearly with the size of the system. As a result, the simulation 
of quantum systems is more efficient on a quantum computer (see, e.g., [lido)) than on a 
classical computer. 

In general, quantum systems are never perfectly isolated from their environments, and 
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in some cases they must be treated as open systems. Understanding the dynamics o: 
quantum systems is therefore important for studying various quantum phenomena [9 
However, the simulation of open quantum systems on a classical computer is also a hard 
task. In addition to the exponential increase in the size of the Hilbert space of the open 
system, one also has to consider the effect of the environment, which adds even more degrees 
of freedom to the problem. 

The Hamiltonian for an open quantum system coupled to an environment can be ex- 
pressed as 

H = H S + H B + H I} (1) 

where H$, Hb and Hi represent the Hamiltonians of the open system, the environment and 
the interaction between the open system and the environment, respectively. The interaction 
Hamiltonian Hi can usually be written in the form 

Hj = Y^Ai® B h (2) 

i 

where the operators A% and Bi act in the state space of the open system and the environment, 
respectively. In general, the environment has a large number of degrees of freedom. Therefore 
the Hamiltonian in Eq. (1) for the global system becomes too complicated to be solvable. 
However, one is usually only interested in the evolution of the open system, and not the 
environment. It turns out that for environments that have short correlation times (i.e., no 
memory effects) and induce slow decoherence in the system, the microscopic details of the 



environment are not important. For purposes of analyzing the dynamics of the system, one 
only needs to know a quantity called the spectral density of the environment. The spectral 
density characterizes the frequency distribution of the noise from the environment on the 
open system. It combines the density of the environment modes and the strength of the 
coupling between the environment modes and the system. For a small set of discrete modes, 
the spectral density would consist of a sequence of peaks. However, the frequencies of the 
environment modes are usually so dense that the spectral density is a smooth function of 
frequency [12 1. 

Then, the task for the study of the dynamics of open systems can be formulated as 
follows: Given the Hamiltonian of the open system, the operators through which the open 
system interacts with the environment, the spectral density of the environment and the 
temperature, how can we obtain the dynamics of the open system? And how can we obtain 
the evolution of the various physical properties of the open system? 

In the case where the environment has no memory effects (also referred to as Markovian 
dynamics), the evolution of the system can be obtained by solving the so-called master 
equation, which in the interaction picture has the form [k^ 



where p$(t) represents the density matrix of the open system, and Ai(u) is the decom- 
position of the operator A± into eigen-projectors of the Hamiltonian The operator 
represents the z-th operator that acts in the state space of the open system in the interaction 
Hamiltonian as shown in Eq. (2). The operators Ai(uS) can be very complicated in matrix 
form depending on the details of the open system. The coefficients Sij (to) are given [l3| by 



Sij (w) 
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where 
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Bi (t) = exp (iHst) Bi exp (—%Hb€) 
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and we have set h = 1. The operator Bi represents the z-th operator that acts in the state 
space of the environment in the interaction Hamiltonian as shown in Eq. (2). The non- 
negative quantities 7^ (to) play the role of decay rates for different decay channels of the 



open system and are given in terms of certain correlation functions of the environment (l3[ 
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dte^{B\{t)BM)- (6) 
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If one were to try and simulate the dynamics of a large open quantum system on a 
classical computer, one would be faced with the problem that the number of basis states 
grows exponentially with the size of the open system. The master equation that describes the 
dynamics of the open system becomes too complicated to be exactly solvable, and sometimes 
it is even practically impossible to derive the master equation. One natural possibility for 
tackling such problems is therefore to use quantum simulation. 

There has been some work on the quantum simulation of open systems in the literature. 
In Ref. 15), the authors suggested an approach for simulating the Markovian dynamics of 



an open quantum system on a quantum computer. They showed that the simulation of 
the Markovian dynamics can be reduced to building generators for a Markovian semigroup. 
However, as mentioned above, the generator for the Markovian semigroup can be difficult 



to obtain for a large system. In Ref. lp |. an approach for preparing the thermal equilib- 
rium state of an open quantum system was suggested. Small-scale open system quantum 
simulators have also has been demonstrated experimentally [n], Q] 

In this paper, by extending the approach of Ref. 16j . we present a quantum algorithm 
for simulating the Markovian dynamics of an open system given the following information: 
the Hamiltonian of the open system, the operators through which the open system interacts 
with the environment, the spectral density of the environment and the temperature. This 
information forms the input of the simulation. 

The structure of this work is as follows: In Sec. [TTJ we present an algorithm for simulating 
the dynamics of an open quantum system. In Sec. IHIl we estimate the resource requirements 
for the algorithm. In Sec. IIVt we provide an example for the algorithm. We close with a 
conclusion section. 



II. THEORETICAL DESCRIPTION OF THE ALGORITHM 



In general, there are three steps in simulating the dynamics of an open system on a 
quantum computer: first, preparing the initial state (see, e.g., |19l-|24j]) of the open system 



and the environment; second, implementing the dynamics on the open system and finally 
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reading out the state of the open system. We will focus on the step of implementing the 
dynamics of the open system, and briefly describe the other two steps of the algorithm in 
subsection II. C since they have been discussed in the literature. 



A. Constructing a model Hamiltonian for the global system 

In general, the details related to the environment in Eq. (1) are unknown. However, 
as mentioned above, one does not need to know all these details in order to obtain the 
system dynamics. One therefore has a good amount of freedom in constructing a model 
Hamiltonian for the environment. One could therefore say that under the Born-Markov 
approximation the spectral density is the only piece of information that one needs to know 
about the environment (see e.g. Ref. [12I, (l^)- Therefore, as long as the spectral density of 
the model environment matches that of the real environment, the effect on the system will 
be the same. 

In theoretical studies it is common to model the environment by a bath of harmonic 
oscillators. For an open system in such an environment, the environment Hamiltonian takes 
the form 

= (*!&* + 5), (7) 

k ^ ' 

and the interaction Hamiltonian can be expressed as 

H I = y ^c k A k ®(bl + b k ^, (8) 

k 

where b\ (b k ) are the creation (annihilation) operators of the environment modes; u k are 
the frequencies of the environment modes and c k are the coupling coefficients between the 
open system and the environment modes; A k are operators that act in the state space of 
the open system and depend on the coupling mechanism between the open system and the 
environment. For a set of discrete modes, the spectral density is usually written as [12I] 

j H = — - w *) > ( 9 ) 

2 ^ m k u k 

k 

where m k is the mass of the fc-th oscillator. The 5 in Eq. (9) is not restricted to infinitely- 
sharp (^-functions, but to 5-function approximants. 

For purposes of simulating the dynamics of an open quantum system on a digital quantum 
computer (which is based on two-state qubits), it is probably more natural to model the 



environment as a bath of two-level systems or spin- 1/2 particles. Such models are also 
sometimes used in theoretical studies (see, e.g., 25|, |26|). For an open system in a spin bath, 



the environment Hamiltonian and the interaction Hamiltonian can be expressed in the form 



B 2 

k 



JZ)"***' (10) 



and 

Ht = 

2 



Hi = \^c k A k ®(g r a x k + g v a z k ), (11) 

k 

respectively, where a k and a z k are the Pauli operators, g r and g^ are coefficients that describe 
the relative size of the transverse coupling and the longitudinal coupling to the open system. 
The transverse component induces relaxation, whereas the longitudinal component induces 
pure dephasing. As will be explained in Sec. II. D, there is a simple alternative method for 
simulating pure dephasing. We therefore take g r = 1 and g v = 0. The spectral density can 
then be expressed as 25] 

J (us) =K^cl8{uj-uj k ). (12) 

k 

The difference between Eq. (9) and Eq. (12) is mostly a matter of convention. 

The environment mode frequencies u k and the coupling coefficients c k can be determined 
as follows: We discretize the frequency spectrum of the environment modes in the full 
frequency range from u min to w max into d elements where each element has a width Aw: 

Ak-'max ^min /i o\ 

CO = . (13) 

Lb 

Correspondingly, the spectral density of the environment, J (u), is discretized and has the 
value J (u k ) for the k-th element, where u k is the frequency of the k-th element that repre- 
sents the k-th environment mode. For a given u k , by making the following approximation 

/ J (cu) du pa J (cu k ) Aoj = nc 2 k , (14) 

JuJk-Au/2 

the corresponding coupling coefficient c k between the k-th element and the open system can 
be obtained. Then the model Hamiltonian is constructed based on the given information of 
the global system and will be used in the implementation of the algorithm. 
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B. Simulating the Markovian dynamics of an open system 



The dynamics of the open system is described by the evolution of the reduced density 
matrix obtained by tracing out the environment degrees of freedom from the density matrix 
of the global system: 

p s (t) = TT B [p(t)], (15) 

where Ps(t) and p(t) are the density matrices of the open system and the global system, 
respectively. The density matrix p(t) undergoes unitary evolution 

p(jt) = U(t,t )p(t )rf(t,t ), (16) 

where 

U{t,t )=exp[-iH{t-t )]. (17) 

Ideally, the dynamics of the open system can be obtained by coupling the open system 
to an environment that has a large number of particles, letting it evolve and reading out 
the state of the open system. In practice, on a digital quantum computer that has a limited 
number of qubits, it is impossible to represent all the particles of a typical environment. 
Therefore we have to use an alternative technique to simulate the dynamics of the open 
system in a large environment. 

The large size of the environment plays a crucial role in justifying the Markovian approx- 
imation. Under this approximation, the typical time during which the internal correlations 
of the environment related to the effects of the open system exist, tr, is much shorter than 
the characteristic relaxation time of the open system, ts', as a result, the influence of the 
open system on the environment is small and can be ignored. Therefore, the state of the 
global system at any time t can be approximately described by the tensor product 

p(t) ^ p s (t) ® p%. (18) 

where is the thermal equilibrium state of the environment and can be written as 

p%=Y,pjm\j = i,-",L=v i , (19) 

3 

with 
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FIG. 1: (a) Quantum circuit for simulating the Markovian dynamics of an open quantum system 
coupled to a time- independent environment. The first register represents the open system and the 
second register represents the environment, (b) Quantum circuit for simulating the Markovian 
dynamics of an open system coupled to a time-dependent environment. 

and \j) (Ej) are the eigenvectors (eigenvalues) of the environment Hamiltonian Hb, Z is the 
partition function 

z = xy**, (2i) 

3=1 

where (3 = l/kgT, is the Boltzmann constant and T is the temperature. Note that in the 
absence of interactions between the environment qubits, the eigenstates \j) can be written 
as a tensor product of the states of the environment qubits 

\j) = |j 2 ) &•••«> \jd), (22) 

and the eigenvalues Ej can be written as a sum of the eigenvalues of the Hamiltonians of 
the environment qubits for the corresponding eigenstate \j) 

E J = E( Jl ) + E(j 2 ) + --- + E( 3d ). (23) 

As a result, the thermal-equilibrium state of the environment is a tensor product of the 
thermal-equilibrium states of the individual qubits: 



Pb = Pt ® Pf ® • - " ® pf, 



(24) 



where 

pt = (l- Pk )\0)(0\+ Pk \l)(l\, (25) 
denotes the thermal equilibrium state of the k-th environment qubit and 

Pk = s—. (26) 

An alternative method that can be used to obtain the Markovian approximation for a 
relatively small environment is to frequently force it back into its thermal-equilibrium state. 
This process can be implemented relatively easily on a quantum computer, where one has 
full access to all the qubits. The procedure for simulating the Markovian dynamics of the 
open system is therefore as follows: Set the environment to its thermal equilibrium state, 
couple the open system to the environment modes and let the global system evolve for some 
time, then reset the state of the environment to its thermal equilibrium state. Repeat this 
process many times and then read out the state of the open system. Mathematically, the 
evolution of the open system in one step is expressed as 

P { t 1] iti + l)r] = Tr B [U(r)fP( 3 r)U\r)\ , j = 0, 1, 2, . . . (27) 

where fP(jr) = p s ® p l B- 

Based on the above analysis, the evolution of the open system can be simulated as 
follows on a quantum computer [see Fig. 1 (a)]: prepare two quantum registers Rg and Rb 
to represent the open system and the environment, respectively; 

(i) on the quantum register Rs, prepare the initial state of the open system; 

(ii) on the quantum register Rb, prepare the thermal equilibrium state of the environ- 
ment; 

(Hi) implement the unitary operation U(r) = exp(-iHr) on the registers Rs and Rb, 
where H is the model Hamiltonian; 

(iv) repeat steps (ii) — (Hi) a number of times. 

(v) read out the state, or the desired observable, of the register Rs. 

In general, the different terms in the model Hamiltonian do not commute. We therefore 
employ the Trotter-Suzuki formula {27] for implementing the unitary operation U(r) = 
exp(—iHr) 
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U (r) = exp [—iHr] 

= lim \ e ~ iH S T / n e - iH B 1 r/n e -iH Il T/n _ _ e -iH Bd T/n e -iH Id T/nl n 

= lim [U 8 (r/n) U Bx (r/n) U h (r/n) .~U Bi (r/n) U Id (r/n)} n 

n—toc 

« [C/ s (r/no) C/ Ba (r/n ) U h (r/n ) • • • ^ (r/n ) C// d (<r/n )P , (28) 

where n is a finite but large number. 

In many cases, one is interested in the value of some physical properties of the open 
system in the thermal-equilibrium state, such as various correlation functions, the partition 
function, etc. The thermal equilibrium state of the open system can be obtained b y rep eating 



the steps (ii) — (Hi) until the open system reaches its thermal equilibrium state 16]. 



Note that the steps followed in implementing the algorithm explained above are the same 



as those used in the algorithm of Ref. 16J. The goal of that work, however, is different 
from the goal of our work. Our algorithm simulates the dynamics whereas the algorithm 
of Ref. [l6[ aims to prepare the thermal equilibrium state of the system. This difference in 
the purpose of the algorithm leads to a number of further differences. For example, in our 
case, we have to design the interaction Hamiltonian and the parameters of the environment 



such that we accurately reproduce the spectral density of the environment. In Ref. [16J, one 
only requires that certain inequalities are satisfied in order for the environment register to 
act as a good environment. In other words, the exact speed of reaching thermal equilibrium 
is not a crucial issue in that work, as long as the thermal-equilibrium state is reached in 
polynomial time. In contrast, if for example there is some symmetry in the simulated system 
that prevents it from reaching the thermal-equilibrium state, then our algorithm would still 
be considered to work successfully if it produces the correct dynamics, even though the 
thermal-equilibrium state is never reached. 



1. Timescales 



At this point we should make a few comments regarding the time scales involved in 
applying the algorithm. The time interval r should be very short compared with rg (so that 
the state of the open system changes slightly during r). The time r can be considered the 
memory time of the environment since the environment is reset to its thermal equilibrium 
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state at every time interval r. Since the timescale for dynamics involving the system and 
a resonant mode of the environment is given by the inverse of q times a matrix element 
between energy eigenstates of the system, the Markovian condition would require that r 
must be small compared to this timescale, such that that the change in the system's density 
matrix is small during the time r. When constructing the model Hamiltonian, the frequency 
spectrum of the environment is discretized into many elements where each element has a 
width Aoj. This width Au must be at most on the order of 1/r to make sure that the 
different 5-peaks have large overlap and produce a smooth spectral density since the width 
of the 5-peaks is on the order of 1/r. 

Note that our algorithm can be straightforwardly generalized to simulate the Markovian 
dynamics of an open system coupled to a time-dependent environment [see Fig. 1 (&)]. In 
each iteration, the state of the environment register is reset to the time-dependent thermal 
equilibrium state of the time-dependent environment. The algorithm can also be generalized 
for the case of an open system that is simultaneously coupled to two or more different 
environments with different temperatures. 

2. Sequential application of the different dissipation channels 

In the above procedure for simulating the dynamics of an open quantum system, we 
have assumed that each mode in the environment is represented by one qubit. In the 
quantum circuit shown in Fig. 1 (a), the environment quantum register R B is prepared in 
the thermal equilibrium state of all the environment qubits with all the different frequencies. 
The operator U{r) acts in the state space of the system qubits and all the environment 
qubits. 

In this subsection, we show that one can reduce the number of qubits required to repre- 
sent the environment by having each qubit represent multiple environment modes. In this 
approach each "evolve- reset" step in the algorithm above is split into multiple evolve-reset 
steps, and in each one of those steps a qubit represents one environment mode. However, 
the qubit is reset to different frequencies in the different steps such that it produces the 
effect of multiple environment modes on the system. Under certain conditions even a single 
qubit can be used to represent the entire environment: the parameters of this qubit are 
sequentially alternated between several different settings such that the sequence covers all 
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the dissipation channels of the environment. 

In the master equation shown in Eq. (3), the derivative of the reduced density matrix of 
the open system is described by a sum over the decay processes of the open system through 
all the different dissipation channels. In the simulation algorithm, the decay of the open 
system in one step of the evolution can be expressed as 

fclti + M = ^B{u int (r)[p s (jr)®p^]U^(r)} 

= (l-rr)pf(jr), (29) 

where J 1 is a superoperator that describes the decay of the open system and it is a sum over 
the contributions from all the different environment modes (the superscript "int" indicates 
the interaction picture). Let 1} denote the superoperator that describes the decay of the 
open system caused by the coupling of the open system to only the j-th environment mode. 
If the open system undergoes a small change during the time r, such that tT <C 1 (i.e., 
all the eigenvalues of r times r are much smaller than 1), then the evolution of the open 
system can be approximated as 

1 _ rF » (1 - rA) (1 - rA) • • • (1 - rr d ) . (30) 

Therefore, the decay of the open system can be simulated by applying the environment 
modes to the open system sequentially. 

As a result, an alternative approach for simulating the dynamics of an open system 
goes as follows: We divide the environment modes into a few sets and let the different 
sets interact with the open system sequentially. In this way, we effectively simulate the 
interaction between the open system and all the different modes in the environment. The 
Hamiltonian for the open system interacting with the i-th set of the environment modes is 
given by 

H® = H 8 + jr i H®+jrH® (31) 

k=l k=l 

where % — 1, 2, • • • , d/di and di is the number of the environment modes in each set. 

One point that requires a little bit extra care here is that each mode in the environment is 
coupled to the system in only one step out of dj di steps, whereas the system Hamiltonian H s 
is applied in all of the steps. One therefore needs to be careful how to calculate the elapsed 
time in the simulated system. If the Hamiltonian in Eq. (31) is applied with a given value of 

12 



r, then after covering all the d/di sets of environment modes the system Hamiltonian would 
have induced a change corresponding to a time r x d/di, while each environment mode would 
have induced a change corresponding to a time r. In order to avoid any problems arising 
from this inconsistency, one can note that decoherence rates are proportional to the spectral 
density (and therefore proportional to c 2 ) from the interaction Hamiltonian. One can then 
use a rescaled Hamiltonian 

H*=H s + j^H% + Mj^H®. (32) 

3=1 V * j=l 

If this Hamiltonian is now applied for a time r, then after covering all the different envi- 
ronment modes both the system Hamiltonian and the coupling to the environment would 
have induced changes that correspond to a time r x d/di, removing the inconsistency in the 
elapsed time. Note that in order to guarantee the validity of the Markovian approximation, 



the rescaled coupling strengths q x y 'd/di must still be small compared to 1/r. 

The new procedure for simulating the Markovian dynamics of open systems is imple- 
mented on a quantum computer as follows: 

i) on the quantum register Rs, prepare the initial state of the open system; 

ii) on the quantum register R B , prepare the thermal equilibrium state of the qubits 
representing a subset of the environment modes; 



-in T 



on the registers Rs and Rb] 



iii) implement the unitary operation Ui(r) = exp 

iv) perform steps ii) — iii) for another set of environment modes, and keep repeating this 
process until the algorithm runs over all the sets of the environment modes; 

v) repeat steps ii) — iv) a number of times; 

vi) read out the state, or the desired observable, of the register R s ; 

In Eq. (30), the error in the decay of the open system introduced by sequentially applying 
the dissipation channels, up to second order in r, is r 2 VV • /i-T,. In the case of using only a 
single qubit to represent the environment, the algorithm will have an error of d 2 T 2 r o , where 
-To denotes the overall scale of the decay rate of the open system caused by the coupling 
of the open system to a single environment mode. When d is large, this approximation 
can introduce a large error. On the other hand, if we use more qubits to represent the 
environment, the errors will be smaller because there will be fewer terms in the sum for the 
error. 
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C. State preparation and readout 



In the algorithm, the initial state of the global system is set to p s (0) <8> p# , where p 5 (0) is 
the initial state of the open system and p 1 ^ 1 represents the thermal equilibrium state of the 
environment. To prepare the initial state of the global system on a quantum computer, we 
prepare two quantum registers Rs and Rb to represent the open system and the environment, 
respectively. In general, the thermal equilibrium state of the environment is a mixed state. 
In order to prepare the mixed state p% as shown in Eq. (19), we can generate a random 
integer j, where j G [1, L], with probability Pj. Then we prepare the corresponding state 
\j) on the quantum register R B . We repeat this step many times. This procedure produces 
an ensemble (in time) of states \j) with the corresponding probabilities Pj. This ensemble 
gives the same effect as the case where the quantum register Rb is prepared in a thermal 
equilibrium state in every time step. 

As discussed in subsection II. A, in the absence of interactions between the environment 
qubits, the thermal-equilibrium state of the environment is a tensor product of the thermal- 
equilibrium states of the individual qubits as shown in Eq. (24). In such cases, the thermal- 
equilibrium state of the environment can be prepared in a simpler way: randomly generating 
or 1 with respective probabilities (1 — pk) and pk, then preparing the corresponding states 
|0) or |1) on the environment qubit, and repeating this step many times, the mixed state 
p^ can be prepared. In this procedure, an ensemble (in time) of states |0) and |1) with the 
corresponding probabilities pk and (1 — Pk) is produced, and it gives the same effect as the 
case where the environment qubit is prepared in a thermal equilibrium state in every time 
step. 



Observing the decoherence dynamics can 
For example, the phase estimation procedure 



)e performed in a number of different ways. 



27| can be used to measure the energy of the 



open system, and one can then monitor the dynamics of the energy distribution. 

Alternatively, any matrix element in the density matrix of the open system can be ob- 



tained using a quantum estimator [28|, |29] as shown in Fig. 2. In the circuit for the quantum 
estimator, if we prepare the second register in the state \ip) and the third register in the state 
\if>), then the value of |(<p|^)| can be estimated by performing single-qubit measurements on 
the index qubit. Through some derivation [28] we have 

^(0) = i(l + |(¥#>| 2 ), (33) 
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FIG. 2: Quantum circuit for a quantum estimator. H represents the Hadamard gate and SWAP 
represents the swap gate. 

where P (0) is the probability for obtaining the state |0) in the index qubit. If the third 
register is in a mixed state p s , then we have 

1 



P(0) 



[1 + (<p\ Ps \<p)) ■ 



(34) 



For the readout of the state of the open system in some chosen basis any diagonal 

element p^ n = (n\p s \n) can be estimated by preparing the second register in the state |n) and 
the third register in the state p s . For the off-diagonal elements p^ n = ( m \Ps\ n )i t ne real part 
of 

Pmn can be estimated by preparing the second register in the state \tp) = (|m) + \n)) / y/2. 
The imaginary part of p^ n can be estimated by setting \ip) = (|m) + i\n)) /a/2. 

The circuit shown in Fig. 2 can also be used for evaluating the expectation values of arbi- 
trary observables. For an operator F, by applying the technique developed in Refs. 28l. |29| . 



one can obtain the value of (F) Ps . Then combined with the quantum circuit for simulating 
the Markovian dynamics of open systems, one can simulate the evolution of the expectation 
values of the physical observables. Various correlation functions can also be obtained using 
this technique. 



D. Simulating pure dephasing of a quantum system 

So far we have concentrated on the case where the nonunitary part of the dynamics is 
caused by the environment modes that are resonant with the open system. This picture is 
valid for energy relaxation. Pure dephasing, on the other hand, is caused by low-frequency 
noise. Such low-frequency noise can also be generated using the algorithm explained in the 
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previous subsections: whenever an environment qubit's state is changed from the ground 
to the excited state or vice versa, the open system feels a telegraph-noise-like change. This 
telegraph noise then causes (mostly) pure dephasing in the system. Although this effect can 
be induced using environment qubits, there is a simpler method to generate telegraph noise. 
The telegraph noise is essentially a classical noise signal affecting the open system. There is 
therefore no need to use qubits in order to produce this classical signal. It can be generated 
using a classical algorithm and added to the system Hamiltonian H$- If different noise signals 
are used in the different runs of the algorithm, the density matrix of the system (averaged 
over the different realizations) will exhibit dephasing dynamics as a function of time. 

For an open system coupled to many ffuctuators, the Hamiltonian for the system can be 
expressed as [30| [see Eq. (11) for comparison] 

H = H s + Y J -Xk(t)A kl (35) 

k 

where 



= (36) 



the random functions ^ (t) characterize the fluctuators' state, instantly switching between 
±1/2 at random times. Therefore, the simulation of the open system coupled to a bath of 
many fluctuators is reduced to simulating a closed quantum system, which will reduce the 
resources required for performing the simulation. 

In simulating the dynamics of a quantum system coupled to many fluctuators, we prepare 
the initial state of the system on a quantum register, then implement the unitary operation 
U = exp (—iHr). Since the signal x {t) is random, to obtain the evolution of the quantum 
system, we have to run the algorithm many times with a different signal x (0 i n each run. 



The noise spectrum of a noise signal is defined as [3J| 



i r°° 

s(u) = - dt coswt<x(t)x(0)>. (37) 

Jo 

In order to generate the telegraph-noise signal [32] , we assume that the environment contains 
a number of fluctuators. Each fluctuator switches between two possible configurations with 
switching rate j i7 and couples to the system with a coupling strength t^. The corresponding 
contribution of a fluctuator to the noise spectrum is a Lorentzian 3l| , 



Si H = - f l ^\ 2y (38) 

47T (u 2 + 7- ) 
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The (low-frequency) noise spectrum of the entire environment is the sum of the noise spectra 
of all the ffuctuators 

By adjusting the parameters r y i and v i: one can produce a variety of noise spectra. A typical 
example in practical situations is 1/f noise. If the number and density of fluctuators is 
sufficiently large, and the distribution of the switching rates of the fluctuators D (7) oc 
and is independent of the distribution of the coupling strengths between the fluctuators and 
the open system, the sum over the fluctuators produces the 1/f noise spectrum (33] 

S(u) = -, (40) 

where G is a constant. 



E. Simulating the non-Markovian dynamics of an open system 

The quantum algorithm presented above for simulating the Markovian dynamics of an 
open system can also be used for simulating a class of non-Markovian dynamics of open 
systems. A common situation where non-Markovian dynamics occurs is the case where a 
small number of degrees of freedom in the environment are coherent enough that they have 



non-negligible memory effects in their interaction with the open system 13j. Examples of 



this situation include the relaxation dynamics of an atom through a cavity that has a high 



quality factor (see e.g., Ref. |l3|) and the relaxation dynamics of a superconducting qubit 
close to resonance with a coherent two-level defect (see e.g. Ref. [3^] ) . 

When a small number of degrees of freedom in the environment are responsible for the 
non-Markovian dynamics, and assuming that one has sufficient understanding of these de- 
grees of freedom (i.e. their intrinsic Hamiltonian, their coupling to the system and their 
decoherence mechanisms and rates), it becomes relatively straightforward to include them 
in the algorithm. One now adds the appropriate number of ancilla qubits needed to describe 
these degrees of freedom. When implementing the evolve-reset part of the algorithm, one 
treats the additional degrees of freedom as part of the system, i.e. they are not reset to their 
initial state. When the measurement is performed at the end of the algorithm, however, 
only the system qubits are used and the additional degrees of freedom are ignored. This 
step corresponds to taking the trace over the state of these degrees of freedom. 
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The introduction of the additional degrees of freedom increases the resource require- 
ments (which will be the subject of Sec. Ill) as follows: the number of qubits used in 
implementing the unitary operation U(r) is increased by \og(D ext ), where D ext is the num- 
ber of degrees of freedom in the environment that are responsible for the non-Markovian 
dynamics. Since the interactions in the expanded system should still be local, the scaling of 
resources will still be polynomial in system+ancilla size and therefore efficient. 



F. Implementing environments other than spin baths 

Our algorithm simulates the effect of the environment using a set of ancilla qubits that 
each represents one spin in a bath of independent spins [351 ] . This type of environment is 
rather straightforward to implement, since the properties and manipulation of the environ- 
ment are done by considering the ancilla qubits one at a time. For a variety of purposes, this 
spin bath is sufficient for the implementation of the desired quantum simulation. However, 
the spin bath has certain limitations. For example, the two-level nature of the environment 
elements (i.e., spins) means that increasing the temperature will reduce the number of spins 
that are in their ground states, thus reducing the relaxation rate induced by this environ- 
ment 26(. This behavior contrasts with the case of a bath of harmonic oscillators, where 
both relaxation and excitation rates increase with increasing temperature. Therefore, if one 
is interested in the temperature dependence of the dynamics, one needs to be careful about 
the differences between different types of environments. 

One possible technique to use a spin bath in order to simulate an oscillator bath and ob- 
tain the correct temperature dependence is to calculate a modified (temperature-dependent) 
spectral density for the spin bath and use this spectral density in the simulation. Alterna- 
tively, different types of environments can be simulated by modifying the algorithm such 
that each element in the environment is encoded into multiple qubits that are treated as a 
single physical object. For example, one could use n environment qubits to represent the 
lowest 2 n energy levels of a harmonic oscillator and then design a Hamiltonian where this 
harmonic oscillator represents one mode in the environment. Thus, one can simulate the 
dynamics of an open quantum system interacting with a bath of harmonic oscillators. Note 
that the state preparation and the form of the Hamiltonian become more complicated in 
this case than in the case of a spin bath. 
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TABLE I: Resource needed for implementing the algorithm. Here [a;] provides the smallest integer 
larger than x, or equal to x if x is an integer. N is the dimension of the Hilbert space of the 
open system, d is the total number of qubits needed in representing the environment in one time, 
di is the number of qubits used in representing the environment in approach 2. m is the number 
of times the unitary operation U(t) = ex.p(—iHr) is implemented, and no is the parameter for 
dividing the time in the Trotter expansion in Eq. (28). 



Algorithm 


# of qubits 


# of operations 


Approach 1 


2[log 2 N]+d+l 


O(mx(2d+l) no ) 


Approach 2 


2[log 2 N]+di+l 


0(mxd/d i x(2d i + l) n °) 



III. RESOURCE ESTIMATION 

In this section, we discuss the resources including the number of qubits and the operations 
needed for implementing the algorithm. 

As shown in Fig. 1 (a), the number of qubits required for simulating the Markovian 
dynamics of the open system is |~log 2 N] + d (here \x] provides the smallest integer larger 
than x, or equal to x if x is an integer), where is the dimension of the Hilbert space of 
the open system and d is the total number of qubits representing the environment. In the 
quantum estimator for the readout of the state of the open system, |~log 2 N~\ + 1 additional 
qubits are needed as shown in Fig. 2. Therefore the total number of qubits for simulating 
the Markovian dynamics of open systems is 2|~log 2 N] + d + 1. One could also use the 
ancilla qubits that are used in representing the environment in the readout of the state 
of the open system, which would reduce the number of qubits used in the algorithm to 
max[[log 2 An +d,2\\og 2 N] + 1]. 

The unitary operation U(r) = exp(— iHr), where H is the model Hamiltonian, is im- 
plemented a finite number m of times. To implement U(r) on the quantum circuit, we 
employ the Trotter-Suzuki formula as shown in Eq. (28), in which U(t) is approximated by 
the product of (2d + l) n ° unitary transformations, where no is the parameter for dividing 
the time in the Trotter expansion. Therefore the number of unitary operations needed in 
the quantum circuit shown in Fig. 1 (a) is m x (2d+ l) n °. Note one could obtain higher 
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efficiency using higher order Suzuki- Trotter formulas, as discussed in Ref. [36] 

In the second approach we presented, the environment is represented using a few or a 
single qubit. The number of qubits required for simulating the dynamics of open systems 
is: 2|~log 2 N~\ + di + 1, where di is the number of qubits representing the environment. 

The unitary operations implemented for each set of qubits are Ui(r) = exp —iH M r , 
where if W is given in Eq. (32). Thus implementing Ui(r) on the circuits by employing 
the Trotter-Suzuki formula requires (2<ij + l) ra ° unitary operations. For d/di sets of the 
environment elements where each set of the elements is implemented m times, the total 
number of unitary operations that need to be implemented is m x d/di x (2<ij + l) n ° . 

In preparing the initial state of the global system, we have to prepare the thermal equi- 
librium state of the environment, which is a mixed state, a number of times. To prepare 
the thermal equilibrium state as shown in Eq. (19), we prepare the state \j) with the 
corresponding probability Pj, and repeat this step many times. In running the algorithm, 
these basis states are fed to the register Rb in the quantum circuit in Fig. 1 (a) one at a 
time, and run the algorithm many times with the the register Rb prepared in the basis 
states. For each basis state \j), we need to reset the state on the register Rb m times in 
the algorithm. To do this, we need to erase the state on the register Rb and then prepare 
Rb in state We can first perform a measurement on the register Rb, the state on Rb 
collapse to a basis state, then we can perform a unitary operation to rotate this state into 
the state 

In the readout of the state of the open system, we employ the quantum estimator, in which 
the information of the open system is obtained by performing single-qubit measurements on 
the index qubit and taking the average. We have to prepare many copies of the state of the 
open system in order to obtain accurate results. Both this procedure and the procedure for 
preparing the thermal equilibrium state of the environment require running the algorithm 
many times. The number of times for running the algorithm, however, does not depend 
on the dimension of the Hilbert space of the open system. Therefore the algorithm can be 
implemented efficiently using 0[m x (2d + l) n °] (or O [m x d/di x {2di + l) no ], if the envi- 
ronment is represented using di qubits) unitary operations. All these results are summarized 
in Table I. 
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IV. EXAMPLE: A TWO-LEVEL SYSTEM IN A SPIN BATH 



In this subsection, we consider the example of simulating the Markovian dynamics of a 
two-level system that is immersed in a thermal bath of independent two-level systems. The 
Hamiltonian of the global system is given by 



k k 



where a x and a z are the Pauli operators, and u s and Uk are the frequencies of the two-level 
system and the environment modes, respectively. The first term is the Hamiltonian of the 
open system, the second term is the Hamiltonian of the environment and the third term 
describes the interaction between the open system and the environment. 

In this ex amp le, assuming Ohmic dissipation, the spectral density of the spin bath is 
expressed as [251 ] 

J (u) = 2irauexp (— u/uj c ) , (42) 

where a is the dissipation coefficient and uj c denotes the cutoff frequency (u c /oj s ^> 1). 
Below we specify frequencies, temperatures and times using a standard unit frequency A . 
We first set uj s /A = 1, a = 2 x 10~ 4 , u c /A = 100, and /3A = 1. The frequency 
spectrum in the region u/A 6 [0.8, 1.15] is discretized into eight elements at frequencies 
co k /A = (0.80 + OMk), with k = 0, 1, • • • 7. The width of each element is Au/A = 0.05. 
The coupling coefficients are determined using Eq. (14). 

For this example, the analytical results for the relaxation rate and the dephasing rate, 
1/Ti and I/T2, can be derived under the Markovian approximation [14] as 

k = \ J{UJ=UJsh k = h (43) 

The Markovian dynamics of the two-level system is simulated with the initial state of the 
open system being set to the excited state |e) of the two- level system. We set the unitary 
evolution time Aqt = 30, and we obtain the evolution of the state of the open system. The 
relaxation dynamics of the two-level system is shown in Fig. 3. 

From Fig. 3, we can see that there is a clear discrepancy between the relaxation rates 
obtained from the numerical simulation and the exact results given by Eq. (43). This 
discrepancy is due to the fact that we used Eq. (14) to determine the coupling coefficients 
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FIG. 3: (Color online) The evolution of the matrix element p ee of a two-level system in a spin bath, 
where p ee denotes the diagonal element of the density matrix that describes the population of the 
excited state. The frequency of the environment mode in the region u/Aq £ [0.8, 1.15] is divided 
into 8 elements with equal width and each element is represented by a qubit. The Hamiltonian for 
the global system is given by Eq. (41). The unitary evolution time Aqt = 30. The red solid line 
represents the analytical result for the evolution of p ee with the initial state |e). The black square 
dots represent the simulated results for the evolution of p ee . 



Cfc even though only eight <5-peaks are used to represent the spectral density: 

k=8 



J (u) = it 2J cl 5 (u - u k ) 



(44) 



Each 5-peak has the analytical form 



k=l 



1 — cos[r (u — uo)] 
7TtA [(u - u )/A ]' 



(45) 
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FIG. 4: (Color online) The relaxation rate of the two- level system as a function of its frequency. 
The black squares represent the results obtained from the simulation of the algorithm; the blue 
dotted line represents the relaxation rate obtained by using the spectral density that is spanned 
by eight <5-peaks in Eq. (44); and the red solid line represents the relaxation rate obtained from 
the analytical spectral density in Eq. (42). Obviously only the first two ones agree well with each 
other. 

In Eq. (14), a good approximation can be achieved when many 5-peaks are used to represent 
the spectral density and the peaks cover a sufficiently wide range of frequencies. The eight 
peaks used in our simulation do not satisfy these conditions. 

In order to further demonstrate the above explanation of the discrepancy, in Fig. 4 we 
show the relaxation rate as a function of the frequency u s of the two-level system. There we 
compare the numerical results obtained from the simulation, the analytical results obtained 
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FIG. 5: (Color online) The relaxation rate of the two-level system as a function of the frequency 
of the two-level system. The black square dots represent the results obtained from the simulation 
of the algorithm using the improved coupling coefficients as discussed in the text; the blue dotted 
line represents the relaxation rate obtained by using the spectral density that is produced by eight 
5-peaks in Eq. (44) using the improved coupling coefficients; and the red solid line represents the 
relaxation rate obtained from the analytical spectral density in Eq. (42). 

from Eq. (43) with the spectral density given by Eqs. (44) and (45) and the coupling coeffi- 
cients determined through Eq. (14); and the exact results obtained from Eq. (43) with the 
spectral density given by Eq. (42). From Fig. 4 we can see that the numerical results for 
the relaxation rate obtained from the simulation fit very well the analytical results of the 
eight-peak spectral density, while both have a systematic deviation from the exact results. 
This systematic deviation from the exact results can be eliminated by plugging in the 
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TABLE II: Results for simulating the Markovian dynamics of a two level system using different 
number of qubits to represent the environment to sequentially simulate a total of eight dissipation 
channels. N denotes the number of qubits representing the environment modes. The exact result 
for T\ is Tf xact = 2/J(u = Aq), and T| xact = 2T±. The results here are obtained using the improved 
choice of the coupling coefficients as discussed in the text. The results are essentially independent 
of N. 



N 


1 


2 


4 


8 


rj-i^ y^exact 


0.998 


0.998 


0.998 


0.996 




1.994 


1.990 


1.991 


1.991 



analytical form of the 5-peaks as shown in Eq. (45) into Eq. (44) for the eight elements 
and then obtaining an improved approximation for the coupling coefficients c&. In Fig. 5, 
we show the relaxation rate as a function of the frequency of the two-level system using 
the improved choice of the coupling coefficients. One can see that in the region around 
u s /A Q = 1.0, the numerical results are now in good agreement with the exact results. 

We also perform the same simulation with the sequential application of the different 
dissipation channels. In different simulations we use different numbers of qubits to represent 
the environment. We use 1, 2, 4, or 8 qubits to represent the environment using the improved 
choice of C/». The ratios Ti/T^ xa,ct and T 2 /T 1 exact are shown in Table II. One can see that 
they are in good agreement with the exact results and that the different simulations give 
essentially the same results. 

We do not perform any numerical calculations for the simulation of pure dephasing using 
telegraph noise here, because such calculations would follow closely similar calculations 
that have been performed in the literature in theoretical studies of telegraph-noise-induced 



dephasing (see, e.g., Ref. [3ll. l37j|). 



V. CONCLUSION 



In this paper, we have presented an algorithm for simulating the Markovian dynamics 
of an open quantum system. The algorithm takes as an input the Hamiltonian of the open 
system, the operators through which the open system interacts with the environment, the 
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spectral density of the environment and temperature. One therefore does not explicitly 
deal with the master equation describing the dynamics. In the simulation, the environment 
is represented by a set of ancilla qubits that are designed to have the same effect on the 
open system as the simulated environment. We have also shown that different dissipation 
channels can be implemented sequentially, thus reducing the number of qubits needed to 
represent the environment. Pure dephasing also allows a reduction in the number of needed 
qubits, since it can be induced by a properly designed classical noise signal. The algorithm 
can also be used to simulate non-Markovian dynamics. 

In the present algorithm, the ancilla qubits play a rather passive role in the sense that they 
only facilitate the dissipative dynamics of the system. These ancilla qubits could be used in 
a more active role as probes or actuators for the open quantum system. By monitoring the 
response of these ancilla qubits as they interact with the open system, one could obtain the 
energy spectrum of the system. Once the spectrum is known, the ancilla qubits can also be 
used to provide or absorb any given amount of energy and guide the system to any desired 
energy eigenstate. The details of this algorithm will be presented elsewhere. 
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